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Abstract 

We consider a class of models describing the dynamics of N Boolean variables, where 
the time evolution of each depends on the values of K of the other variables. Previous 
work has considered models with dissipative dynamics. Here we consider time-reversible 
models, which necessarily have the property that every possible point in the state-space is 
an element of one and only one cycle. 

As in the dissipative case, when K is large, typical orbit lengths grow exponentially with 
A, whereas for small enough A, typical orbit lengths grow much more slowly with A. The 
numerical data are consistent with the existence of a phase transition at which the average 
orbit length grows as a power of A at a value of K between 1.4 and 1.7. However, in the 
reversible models the interplay between the discrete symmetry and quenched randomness 
can lead to enormous fluctuations of orbit lengths and other interesting features that are 
unique to the reversible case. 

The orbits can be classified by their behavior under time reversal. The orbits that trans- 
form into themselves under time reversal have properties quite different from those that 
do not; in particular, a significant fraction of latter-type orbits have lengths enormously 
longer than orbits that are time reversal-symmetric. For large A and moderate A , the 
vast majority of points in the state-space are on one of the time reversal singlet orbits, and 
a random hopping model gives an accurate description of orbit lengths. However, for any 
finite A, the random hopping approximation fails qualitatively when A is large enough 
(A » 2(2^)). 

Key words: Gene Regulatory networks; Random boolean networks; Time-reversible 
Boolean networks; Cellular automata; 



1 Introduction 



1.1 Review of dynamical boolean networks 



In recent years considerable effort has been devoted to the study of the development of 
complexity in dynamical systems. Complexity is observed in such different examples 
as ecosystems (see [1]), spin glasses (see [2]), and quite broadly through the biological 
sciences [3]. One thread of activity involves the study of the behavior of dynamical 
systems consisting of variables al (the site label j = 1, . . . , A^), where each o"^ 
is Boolean, so that it takes on one of two values that we choose to be ±1. The 
configuration of the system at time t is characterized by a 'state' E^: 

E, = {al,al...,ai,...,an. (1) 

The time development is given by saying that the state at time t + 1 is a prescribed 
function of the state at time t, i.e. 



Et+i = M{^t) ■ (2) 
The mapping can also be written 

al, = F^{Et) for j = l,...,N (3) 

where the also take on the values ±1. 

The number of different states of the system is finite; it is 

c^ = 2^, (4) 

and the mapping function Ai is independent of time. Therefore, starting from any 
state, eventually the system falls into a cyclic behavior and follows that cycle forever. 

Typically, each of the functions F-^ (E) is picked so that it depends upon exactly K 
distinct input spin variables in the vector E. A random choice is made to determine 
which components, will appear in each F-'(E); this fixes the 'wiring' of the real- 
ization. Once K and A^ and the assignment cr'^'s are fixed, the each mapping function 
F^{T,) is selected at random from the set of all Boolean functions of K Boolean vari- 
ables. The randomly chosen N sets of K input spin variables and the functions F^ 
compose a realization. Given an initial configuration of the system, a realization will 
completely define the system's behavior. The realization is picked at the beginning 
of the calculation for the system and remains independent of time. Since there are 
2^^ Boolean functions of K Boolean variables and different ways to assign K 
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distinct input spin variables, as K and/or N become large the number of different 
realizations is truly huge. 

This kind of model is often called a Kauffman net because Stuart Kauffman[3-5] 
developed a program of study for generic maps of this type. Later, this program 
was extended by Dcrrida[2], Flyvberg[6], Parisi[7-9], and others[10,ll]. Quantities 
of interest that have been studied include the distribution of cycle-lengths and the 
number of starting points which will eventually lead to a given cycle. (These points 
form what is called the basin of attraction of the cycle.) One calculates the above 
quantities by enumeration or by Monte Carlo simulation for each realization and then 
averages over realizations with the same values of N and K. 

The Kauffman net is said to be 'dissipative' since several different states may map 
into one. Thereby information is lost. 

The behavior of Kauffman nets arc interesting and surprising. Large values produce 
a complex time-behavior which closely resembles Parisi's theory of spin glasses [2, 12,8]. 
For large K, the cycle lengths grow exponentially with the number of spins A^". Con- 
versely, for K = OT K = 1, the cycles tend to be short. At the 'critical' value, 
K = 2, typical cycle lengths grow as a power of A^, and for large the probability of 
observing a cycle of length L varies as a power of L. (For a study of critical properties 
see Refs [8], [10] and [13].) This three-phase structure is typical of phase transition 
problems[14] in which an ordered and a disordered phase are separated by a critical 
phase line. 

1.2 A time-reversible network 

Thus, a great deal is known about the behavior of Kauffman nets, which can be viewed 
as a class of generic dynamical mapping problems. But not all problems are generic. 
For example, many of the systems considered in Hamiltonian mechanics are reversible. 
Such systems have the property that some transformation of the coordinates (for 
example changing the sign of all velocities) makes the system retrace its previous path. 
Thus a forward motion and its inverse are equally possible. This paper is devoted to 
a study of the behavior of discrete reversible maps. 

In contrast to a dissipative system, in a finite and reversible dynamical system every 
possible state is in exactly one cycle. Because one and only one state at time t maps 
into a predefined state at time t + 1, any cycle can be traversed equally well forward 
or backward. There are no basins of attraction in this kind of system. The long-term 
properties are then described by giving the number of cycles of length I, N{1). 

It turns out that for the smaller values of K, N{1) is an oscillating function of / in 
which the small prime divisors of / play a major role. We shall study this effect in a 
companion paper. For now, we focus on the gross scaling properties of N{1), by using 
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a cumulative distribution 

which gives the probabihty of finding a cycle of length greater than / by picking 
the realizations and the initial cycle element at random. We call S{1) a "Survival 
probability" . 



1.2.1 Definition of the model 

We construct our time-reversible maps using the method of two time-slices that was 
studied for discrete systems by Fredkin and collaborators.Q[15] 

As in the usual dissipative Kauffman model, our basic variable is a list, S, of the N 
'spin variables', a^: 

S = (a\a2,...,a^...,0. (6) 



The value of S is given for each integer value of the time, and written as S^. In our 
reversible mappings the spin configuration at time t + 1, T^t+i, depends upon and 
Si„i according to the rule 

^Ui = ^i-iF^i^t) , (7) 

where the F's are picked exactly as in the dissipative Kauffman net. Since the cr's 
take on the values ±1, our model can be written in the equivalent form 

al,al, = , (8) 

which exhibits a quite manifest time reversal invariance. The models given by equa- 
tion (7) are the subject of this paper. 

The information needed to predict future time steps is called the state of the system. 
In our case, the state at time t, St, is given by two time slices or substates and 
as 



^ The construction of time-reversible models by using variables from two different times has 
a long history. Imagine doing a calculation in ordinary classical mechanics using small but 
discrete time-steps. The behavior depends upon the position and velocity of each particle, 
but one might wish to do the analysis in terms of positions alone. To do this, one works with 
a state defined by the positions at two closely neighboring times. The difference in position 
at the two times gives an estimate of the velocity vector. In this way one can construct a 
two slice model of particle behavior. 
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St 



) 



(9) 



and the full mapping is of the form 



(10) 



The history of the system is given by listing the substates in order as 



S: 



'2: 



'n 



If there are N spins, the volume of the state space is 



= 2^^ = 0;^ 



(11) 



Since the dynamics are deterministic, fl gives the length of the longest possible cycle. 
However, as we shall see, that length is never attained. 

1.3 Questions to be asked 

This paper concerns the distribution of cycle lengths in the time-reversible models. 
Figure 1 shows plots of probabilities of observing a cycle of length larger than / for 
systems with N — 10 and various K, averaged over reahzations. For K — 1 the cycles 
are very short; for K = N they have a wide range of lengths, but the longest ones 
have length of order 2^, much less than the number of points in the state space, 
2^^. For K = 2, a wide range of cycle lengths is seen, inchiding some lengths which 
considerably exceed the ones in the K = N system. For smaller values of K, the 
number of cycles of length I tends to be a strongly oscillatory function of I. We shall 
discuss this oscillation in a subsequent publication. 

To understand these results for the cycle lengths, we shall need to understand in some 
detail the interplay between the quenched randomness of the system and the time- 
reversal symmetry. We will find that the cycles can be divided into two classes; those 
that are symmetric under time-reversal and those that are not. There are profound 
differences between the behaviors of these two types of cycles. 

The companion paper will discuss the growth of the 'Hamming distance' between 
configurations, defined as the the number of Boolean variables which are unequal in 

the two configurations. We shall compare two configurations that initially differ by 
a single spin fiip, and see how this distance depends upon the number of iterations. 
For small K, the growth in the distance may be very slow; for large K we see much 
more rapid growth. 
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Fig. 1. Plots of probabilities of observing a cycle of length greater than I as a function of I 
for A'' = 10 and K = 1,2, and N. These plots are derived from simulations which average 
over 10, 000 realizations for each value of K. For each realization, one initial state randomly 
chosen from the state-space was examined. 

Our simulational data of cycle length distributions suggest but do not prove the 
existence of phase transition. However, in the companion publication wc shall see 
that the behavior of the Hamming distance can be used to demonstrate convincingly 
the presence of a type of percolation transition at a value of K of about 1.6. 



1.4 Outline of the rest of paper 



In the next section, we discuss some special features of time-reversible models and 
their implications for classifying different kinds of cycles. The section after that is 
devoted to the limiting cases, K = and N. Section four describes the structures 
seen at intermediate K. The appendices cover some more peripheral issues. 
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Mirror 



Fig. 2. Mirrors produce a closed cycle. The mirrors are at S3 and between the two Si's. 
Notice that interchanging the direction of the arrows gives the same result once again. 

2 Time reversal invariance 



At first sight, it is not obvious that time- reversal invariance should have any important 
effect on the distril:)ution of cycle lengths. However, G. Birkhoff[f6] and later John 
M. Greene[17] and others[18] showed one important mechanism by which the time- 
reversal process works to set the cycle length. 



2. 1 Symmetry points 



In our model there are special points in the phase space which we might call mirrors. 
A mirror produces a time-reflected motion in the sequence, for example 

. . . , E3, E2, Si, El, E2, S3, . . . (12a) 
or, as another example 



. . . , E3, E2, El, So, El, E2, E3, . . . (12b) 

We call the first type (equation (12a)) a twin configuration and the second type 
(equation (12b)) a sandwzc/i configuration. A sketch of the mirrors is given in figure 2. 
The phase space contains many of these mirrors. A typical and important cyclic 
motion is for the sequence of substates to hit a mirror, be reflected, hit another 
mirror, be reflected once more, and thereby be forced into a cyclic behavior. 

More explicitly, if we record an orbit of the time-reversible model using a sequence of 
substates, for example . . . , Et_i, E^, E(_|_i, . . . , we may call the state 5- value defined 
as (s^' J at which S^+i = E^ a twin special point, and another 5-value a sandwich 
special point, if St+2 in the sequence equals S^. Thus a twin point is any state of the 
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form 




so that the cycle will spread out in a palindromic fashion before and after this special 
point in the pattern of equation (12a). Since the entire volume of the state space is 
r2 = a;^ and since there are uj of these invariant points, the chance that a randomly 
chosen point in the state space is a twin point is 1/uj. Correspondingly, a sandwich 
point appears when there is some value of T,^^^ for which 

= 1 for all j (14) 

(see equation (7)). A sandwich point gives rise to a sequence of the form of equation 
(12b). The number of sandwich points depends on the realization. For example, if 
= —1 for any j, then there are no sandwich points at all. For a given realization, 
we denote the number of substates that have the property of equation (14) by m, and 
we denote each such substate by the symbol E^^^; the a label differentiates between 
the m different sandwich substates. Since for any substate F a state of the form 
S — (j,o>)) is a sandwich state, there must be mu of these sandwich states. 

Appendix A discusses the properties of the sandwich points. As shown there, the 
average over realizations (m) is unity. When K — N and N is large, m follows a 
Poisson distribution, 

p(m)^-^ (N^oo). (15) 

m\e 

When K is small, there are large realization-to-realization fluctuations in m, and the 
average over realizations of higher powers of m, for example (m^), can grow rapidly 
as a function of N. In fact, when » 2^^^^ the vast majority of reahzations have 
m = 0. 



2.2 Inversions and cycles 

Two states S and S' are time reversed images of one another US— while 

S' — {^^- All cycles belong in one of two classes. The first class, which we call 
special cycles., each contains at least one pair of time-reversed images of one another. 

The other class, called regular cycles^ contain no such pairs. In Appendix B we show 
that each special cycle of length greater than one contains exactly two distinct special 
points. Furthermore, if the points are both twins or both sandwiches, the cycle has 
even length; if they are of different types, the cycle length is odd. 
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2.3 Counting special points and special cycles 



The arguments in subsection 2.1 imply that for a given reahzation there are mw 
sandwich points and uo twin points. Among them there are m points that are both 
twin and sandwich points (consisting of a sandwich substatc followed by the same 
substate), all of which produce m cycles of length one. Thus, there are (m + l)uj — m 
distinct special points in the state space. The total number of special cycles is the 
number derived from the points which are both twins and sandwiches, m, plus the 
number derived from all the other special points, [(m + l)u! — 2m]/2, yielding 



number of special cycles = - — ^ . (16) 

The importance of these special points in determining cycle properties will become 
clear in the following sections. 
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3 Limiting cases 



This section discusses the behavior of reversible Boolean nets for the cases K — 
and K^N. 

3.1 K^O 

When K = the evolutions of the different cr-^'s are uncorrelated. Each aj repeats 
after either one, two, or four steps. (In contrast, in the Kauffman net after the first 
step each of the cr-^ 's remains constant.) For large N, each system is likely to contain 
at least one with period four. Hence cycles of period four will dominate. There 
will be roughly N{4:) a;^/4 such cycles. 

The uncorrelated cycles of the spins produce a hamming distance which can have 
value zero or one and has period one, two or four. 

3.2 K^N 

In this part we first review the results for the dissipative Kauffman net 
and then proceed to discuss the time-reversible case 

3.2.1 Dissipative case 

The case in which K has its maximum value, K—N, was first analyzed for the dis- 
sipative Kauffman model by Derrida[2]. This case has the simphfying feature that a 
change of a single spin changes the input of every function FK Since the functions 
are chosen randomly, this means that every new input configuration leads to an 
output configuration Ej+i that is picked at random from the whole phase space, with 
its volume 

a; = 2^. 

This process continues until the time T at which = E^ for some r < T. After 
that, the system cycles repeatedly through the sequence E,-, . . . , E^-i. 

To find the distribution of orbit lengths, we first calculate Pn, the probability that 
starting from a randomly chosen initial state at time t — 0, the orbit closes at time 
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n. To do this, we define to be the probabihty that the cycle remains unclosed after 
n steps. The probability that a closure event occurs at time zero is po = 0, and thus 
go = 1- At the time one, the system has a probability pi = l/u of falling into the 
initial value and a probability gi = 1 — pi of not doing so. At time two, there are two 
possible cycle-closures, since the new element can be the same as either the zeroth or 
the first element. Thus the conditional probability of a closure at time two, given that 
the closure event did not occur at any earlier time, is 2/lj. Similarly, the conditional 
probability of a closure event at time t = n, given that the system has not closed at 
time t = n — 1, is just n/uj. Thus the likelihood of a cycle closure at step n is 



Pn = -Qn-i (17a) 



and correspondingly the g„ satisfy 



n 



qn= 1 Qn-l ■ (17b) 



The solution to equation (17b) with go = 1 is 



9n=n(l-i'M- (18) 

We shall see that g„ -C 1 unless n uj. Therefore, we can write 



lng„ = X:in(l-jM (19) 
and expand the right hand side for small j/uj, yielding[^ 

qn = exp[-n{n + l)/2uj] . (20) 
The probability Pn of obtaining a cycle closure at time t = n is then 



Pn = -qn^i = ^e-"("+y/2- . (21) 

UJ LJ 



To obtain V{L), the probability that a given starting point is in the basin of attraction 
of a cycle of length L, we note that a closure event at time t = n yields with equal 
probability all cycle lengths up to n. Therefore, 



^ One can also write g„ = uj ^lj\/{uj — n)l and expand the factorials using Sterling's 
approximation. 
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oo 



Pn 



n 



(22) 



n=L 



which is well-approximated by 



I" 

V{L) ^ / 




1 + 2L 



) 



) 



(23) 



(24) 



Note that the probability distribution of cycle lengths is asymptotically Gaussian (for 
large lengths), and that the length of a typical cycle is of order ^^[5,8]. Derrida has 
pointed out [19] that the random hopping assumption (or annealed approximation) 
is exact for the dissipative Kauffman net with K = N . When K < N, the random 
hopping assumption is no longer exact. However, the K — N results agree well with 
the simulational data for large but finite K. 



3.2.2 K ^ N: Reversible case 

As we shall see, the behavior of the reversible model in certain parameter regimes 
is also well-described by a model in which the time-development can be considered 
to be random until a closure event occurs. However, the analysis oi K — N limit is 
more subtle for the reversible model than for the dissipative case. 

Wrong calculation: leave out special points. To illustrate some of the complica- 
tions that arise when we consider the reversible model, we first present a naive (and 
wrong) adaptation to the reversible system of the argument in subsection 3.2.1. Note 
that in the reversible models, the state at time t, St, depends on the spin configu- 
rations, or substates, at two times, S^.i and S^. We once again consider a sequence 
of states SojSijS2, ■ ■ ■ and assume that the map induces "random hopping" through 
the state space (each Sj chosen with equal probability from all allowed possibilities). 
If the state <S„ happens to be the same as Sq, then the cycle closes, with the cycle- 
length being L — n. Note that at the nth step there is only one possible output that 
will give closure, Sn = Sq\ the n — 1 other values of Sj, 1 < j < n, arc impossible 
because each cycle must be traversable both forward and backward. Therefore, if the 
cycle has not closed in the first n — 1 steps, the total number of allowed possibilities 
for Sn is — n, of which only one will give closure. This argument yields an estimate 
for Pn, the probability of closure at the nth step, given that the orbit has not closed 
previously: 



1 



(25) 



Pn = 



— n 
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Therefore, this estimate imphes that pi, the probabihty that the orbit closes at the 
Ith step, should be 



in the limit of large Q. 

However, equation (27) is wrong. Looking back at figure 1, one sees that for K = 
N — 10, the average cycle length is of order a; = 2^ pa 10"^. However, equation (27) 
implies an average cycle length of order = 2^^ 10^. 

A more accurate calculation. The problem with equation (27) is that we have 
ignored the role of the special points. A sequence of substates of the form 



is reflected at the twin point and must continue Sfj+j,. . . S^j. Similarly, a 

sequence of substates of the form 



is reflected at the sandwich point. After the first special point has been hit, the orbit 
retraces and then continues until a second special point is reached. Once the second 
special point is reached, say at time t = t*, the orbit is reflected again, and closure in 
less than t* additional steps is guaranteed. Since twin points are hit with probability 
l/uj and sandwiches with probability m/u at each time step, this mechanism yields 
orbit lengths of order u; rather than the 0{fl) result of equation (27). 

In the dissipative case discussed above, when K = N, the calculation of cycle lengths 
is exact. We have not been able to do that well in the reversible case. However, we 
present here a simple approximation for the distribution of cycle lengths that is rather 
accurate when cu is large and I is less than or of the same order as cu. 

We begin by specifying a realization of the network. Given this realization, we consider 
a sequence of substates 

G/ = Eo, El, E2, Ei, E;+i (28) 
produced by the map. 




(26) 



k=l 



exp[— Z/fi] 



(27) 
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We define a regular sequence to be one which has neither closed nor reached a special 
point. A regular sequence can be used to construct a part of either a regular or a 
special cycle. We define the sequence (28) to be a sequence of length /. Given a regular 
sequence Gi, we may produce a by evolving Gi for another step. We define the 
probability q{l) as the fraction of the realizations in which Gi+i is also a regular 
sequence and the probability p{l) as the fraction in which it is not. 

We take Gi as a regular sequence and now imagine calculating the next substate 
E;_|_2. As an approximation, assume that all E;+2 appear with equal probability, -. 
The probability that E;+2 — "^i+i (i-e., ^i+i is a twin point) is 1/cj, and the probability 
that E;+2 = S/ (i.e., is a sandwich point) ism/u. There is also a chance l/u that 
E;+2 = Eq. In this last case, the orbit will close with no special points if, in addition, 
E;+3 = El. Thus, this estimate yields a probability of closure without special points 
that is of order as in our naive estimate above. 

Therefore, p(/), the probability of a closure event at step /, given that the sequence has 
not closed previously, is the sum of two terms: the probability for closure to a regular 
sequence ~ 1/0, and the probability of getting a special point p*^ ~ (m + l)/a;, 
so that 



In addition to ignoring the possibility that Ei+2 has already appeared in the sequence, 
we have ignored terms of relative order l/u. This is quite reasonable for large N. 

Now we can derive expressions for the number of cycles of different types. There are 
Q different starting configurations for sequences. We require that the first substate 
Eq not be a sandwich substate and that the second substate Ei not be equal to Eq, 
Therefore, the fraction of sequences of the form Eq, Ei (e.g., / = 0) that are regular 
is (1 — l/uj){l — m/uj) ~ (1 — (m + l)/tt;). Each iteration reduces the fraction of 
sequences that are regular by a factor of 1 — p, until after / steps we find that the 
number of regular sequences, Njis{l), is (where again we disregard terms of relative 
order I/lu and smaller) 



Since the probability of closure to a regular cycle is l/fl, the probability that a 
randomly chosen point in the phase space is part of a regular cycle that closes in I 
steps, Pr{1), is Pr{1) = Njis{l)/^. Because each cycle of length I is found by starting 
at any of / points on the cycle, the average number of regular cycles which close after 
/ steps is 



p{l)^ p'^{l)+p^{l) « {m + l)/uj . 



(29) 



Nrs{1) = J^(l - (m + l)/ou){l - pY ^ fle-^'"+^^^/' 



(30) 



Nr{1) 



Nrs{1) 
ni 



^-lg-(m+l)//w 



(31) 
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A very small proportion of the points in the state space are, in fact, parts of regular 
cycles. The number which take part in regular cycles of all lengths, Mr, is 



y^R = T.jMj)^-—T- (32) 

This number is indeed much smaller than the state space volume 

We now turn our attention to the special cycles, which dominate the state space in 
this high K limit. All of them can be found by starting at a special point and iterating 
until a second special point is reached after I steps. Then the orbit reverses itself and 
closes after an additional / steps. There are three kinds of special cycles, twin-twin, 
sandwich-sandwich, and twin-sandwich. There are uim different ways to choose the 
initial point if it is a sandwich point and uj ways to choose it if it is a twin point, so 
the number of twin-twin cycles and sandwich-sandwich cycles of length / are: 

Nu{l) = le-('"+i)^/(2c.) (33^) 

iV,,(/) ^ !^e-('"+^)'/(2'^) . (33b) 

The factors of two arise in equation (33) because each cycle of these types is found 
twice by this method. Similarly, the number of odd (sandwich-twin) special cycles is 



Note that this estimate for the distribution of special cycle lengths depends exponen- 
tially on Z, in contrast to the Gaussian dependence for the dissipative model. 

To check these conclusions we plot in figure 3 the number of cycles of length 

/ as a function of I averaged over realizations. The realizations used all had m = 0. 
The theoretical estimates (equations 31 and 33) agree very well with the numerical 
results. 

Average over m. Our results of equations (31) and (33) depend upon the number 

of sandwich special points, m. In this section, we shall denote averages over m by 
< ■ >. In Appendix A we show that for K — N , the probability distribution for m is 
a Poisson distribution 



pirn) = ;— with A = 1 . (34) 

m! 

Averaging equation (31) using the weight defined in equation (34) gives that the 
reahzation average of the number of regular cycles of length I, < Nr{1) >, is 
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Fig. 3. iV(/) as a function of I, where N(l) is the number of cycles of length I, averaged over 
realizations. The curves are the results of the random hopping approximation (equations 
(31) and (33)) and the points are simulational data. Simulation results are plotted after 
being averaged over 10,000 realizations with K = N = 10 and m = 0. 



< Nr{1) >= r ^ exp [(e"'/'^ - 1) - l/uj 
Similarly, for the various kinds of special cycles 



(35a) 



< Nu{l) > = ^ exp [(e-'/(^'^) - 1) - //(2a;) 

< Nts{l) > = exp [(e-'/(2^) - 1) - l/u 



<Ns,{l)> 



- exp 



"(6-^(2^^) -1)-// 



(35b) 
(35c) 

(35d) 



To test equation (35) against simulations, in figure 4 we plot Ns{l), the number of 
special cycles averaged over realizations, against /. As discussed in subsection (2.2), 
special cycles formed by two twin points or two sandwich points have an even length, 
while the ones with one twin point and one sandwich point have an odd length. 
Therefore, Ns{l) =< iV,,(/) > + < Ntt{l) > if Hs even and Ns{l) =< Nts{l) > 
if / is odd. We find that Ns{l) oscillates because of this difference between even 
cycle lengths and odd cycle lengths, leading to a two- branch structure in Ns{l). This 
structure was indeed observed in our simulations. 

The random hopping analysis of the K = N situation gives a Hamming-distance 
behavior which is exactly the same in the dissipative [5] and reversible systems. In 
both cases, the systems start from a pair of configurations which differ in the value 
of one spin at time t. By the next time-step, the configuration will be random, so 
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Simulation Results 
Predictions of Eq(35) 




even cycles 



10 15 20 25 30 35 

Scaled Cycle Length 1/2 



Fig. 4. The average number of special cycles of length I, Ns{l), plotted against scaled 
cycle length 1/2^ . The result is averaged over 10,000 randomly selected realizations for 
K = N = 10 and hence is an average over m. For each realization all special cycles were 
enumerated. The symbols are the results of the simulations. The theoretical results for 
Ns{l) obtained from equation (35) are included as solid lines. 

that half the spins will be "wrong" . Thus the Hamming distance immediately comes 
to the value N/2, and stays there. 
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4 Intermediate K 



Section 3 outlines both qualitative and quantitive pictures of the limiting cases K — Q 
and N. In this section, we will discuss the system's behavior for intermediate K values. 
As K is increased, there is an evolution from dynamically independent clumps of spins 
to a situation in which there is random hopping over the whole state-space. 

Figure 5 shows numerical results of Ssil) and Sr{1), the survival functions for special 
cycles and regular cycles separately, for various intermediate K values. Recall that a 
survival function of cycle length / is defined to be the probability of observing cycles 
with length greater than I (see equation (5)). Nonintcgcr values of K are produced 
by mixtures of spins with the neighboring integer i^-values. The survival functions 
for special cycles, Ss{l), as shown in part (A), shows no surprising structure. As K 
decreases from 2, the observed cycle lengths get shorter and shorter. For K >3, the 
survival probabilities approach in a uniform manner the result of equation (35). The 
regular cycles are different. For K > 5 wc sec a uniform approach to the K = N 
result. At = 1.6, there seems to be a power law behavior with Sji{l) ~ l^^. Then 
for larger values of K, Sr{1) seems to decrease faster than power law. Finally, K — 3 
and K — A show a remarkable bump at large values of indicating that there is 
a new process going on at large I. 

The difference between part (A) and (B) implies that the relative importance of the 
two pieces of cycle-closing mechanism depends strongly on K. Moreover, the relative 
numbers of regular and special cycles arc also strongly i^-dcpendent. This point is 
illustrated in figure 6, which shows the ratio of the number of regular cycles to that 
of special cycles as a function of K. For small K, several subsets of the spins may 
become uncoupled, and most cycles are regular cycles. Conversely special cycles are 
more likely for large K. 



4-1 Anomalously Long Regular Cycles for Intermediate K Values 

Notice the bump in figure 5 part (B) for X = 3 and K — Am the region 1 <^l <^uj. 
This bump arises from a group of regular cycles which are anomalously long. More 
careful study shows that in fact regular cycles split into two groups, which scale 
differently. To demonstrate this effect, we plot in figure 7 the survival functions Sr{1) 
for K = ?> and for a range of N values against the cycle length / normalized by 
2^. Manifestly the distribution of short cycles scales with u — 2^. This scahng 
can be explained in the following way. There are {m + 1)2^ special points in the 
state-space of the system, as discussed in section 2.3. To obtain a regular cycle, the 
system must not hit any of the special points. The probability of hitting a special 
point ~ 1/2^, so the lengths of the regular cycles scale as 2^. For a more explicit 
and quantitative reasoning, we may appeal to the regular cycle result based on the 
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Fig. 5. The plot of survival functions of scaled cycle length 1 = 1/2 for both special and 
regular cycles for various K values and N = 10. (A) shows the survival function for special 
cycles, Ssil), and (B) shows the survival function for regular cycles, Sji{l). The functions 
were averaged over 25600 realizations for each K value. The survival functions converge to 
the random hopping case for larger K values. However a qualitative theory is still needed 
to predict the function forms for small K. 

random hopping assumption described by equation (35), which provides a reasonably 
accurate approximation of the system when the system is sufficiently well-connected. 
Note that in equation (35), the regular cycle distribution clearly scales as 2^. 

While we have presented an intuitive explanation for the scaling of the relatively short 
regular cycles, one may still wonder why there is a population of anomalously long 
cycles and why it appears and disappears as K increases. To attack these issues, a 
number of realizations which produce extra long cycles were studied and similarities 
among them were observed. In general, a realization and the initial configuration 
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Fig. 6. The ratio of the number of regular cycles to the number of special cycles for N=6, 
8, and 10. Both the numerator and denominator are averaged over 25600 realizations. At 
small K there are many more regular cycles than special cycles; at large K most of the 
cycles are special. 
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Fig. 7. Plot of Sji{l) against / = 1/2 for K = 3 and various A''. The drop of Sr{1) 
corresponding to the relatively short cycles scale as 2^ , which supports the argument that 
they are due to special points. 

have to satisfy the following two conditions to be able to generate an anomalously 
long cycle: The realization should not have any sandwich points, thereby annulling 
the mechanism to close a cycle using them, and the functions assigned to the spins 
together with the choice of initial configurations should prevent a system from hitting 
a twin point. For small K values, one can easily find realizations without any sandwich 
points. For instance, if any spin is assigned the function that is for all inputs, 
then the realization has no sandwich points. In fact, in Appendix A, we prove that 
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almost all realizations have no sandwich points for finite K and large enough N 
(»2(2")). 

The twin points in the cycle could be avoided in many ways. For example, when a 
function assigned to a spin equals to the constant +1, if this spin starts from Uo = +1 
and o"! = —1, it will follow the progression of +1, —1, +1, —1, ... , and the system 
never hits a twin point since at 7^ (Tt+i always holds. Since there is a probability 
1/2*^2 ) lYioX any given spin is assigned the function and a probability 1/2 that 

this spin starts with uq — —ai, the probability p2 that at least one spin is in the 
"H 1 — " progression is 

«^l-exp(-^(^) (2(^")»1). (36) 

Thus, when 3> 2^2^)^ almost all starting configurations and realizations will not 
hit a twin point. 

Another way to avoid twin points consists of two spins that are assigned identical 
inputs and functions. In this case, it is possible to choose an initial configuration for 
these two spins, which will stop the system from forming a twin point. For instance, 
a system in which spin 1 and spin 2 are assigned identical input spins and functions 
starting from cr^^Q = a}^^ — 1 and (t|^q = —1, a^^^ — 1 can never hit a twin point. 
It is also easy to prescribe certain simple progressions for a few spins and stop the 
system from hitting a twin point. For example, two spins both with period 3 evolving 
following the pattern 



(To: + + - + + - + + 

cTi :- + + - + + - + + ■■■ 

can prevent the system from forming a special cycle. 

The mechanisms preventing the system from hitting a twin point always appear 
to be related to some small piece of the system that evolves independently from 
other parts of the system. The couplings are such that this piece is not affected 
by anything outside itself (though in general it can and docs affect the rest of the 
system). We call a piece like this a local structure. Local structures arc to be discussed 
in detail in the next paper in this series. Note that local structures are very unlikely 
for sufficiently large K values, where the spins of the system are quite correlated, 
unless N is enormous. When K is small, local structures occur much more frequently. 
When many a local structures are present, almost all the realizations and initial 
configurations have no special points. Thus the presence of local structures leads to 
very long regular cycles. Since for a given N the local structures become less probable 
as K increases, it is natural to find the number of regular cycles decreases as K. 
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Fig. 8. Survival functions of regular cycle lengths in systems with all functions chosen 
at random as well as systems with one spin assigned the function —1 and a second spin 
assigned the function +1 and initial condition (+1,-1), plotted against for various 

N and for K = 3. As expected, the lengths of anomalously long cycles in the systems with 
randomly functions, as well as all the cycles in the (+, +, — , — ), (+, — , +, — ) systems, scale 
approximately as 2^^. 



Section 3.2.2 presents a naive theory of hopping in state space that ignores the role 
of the special points and predicts that the distribution of orbit lengths in a system 
of N spins should scale as 2^^. This naive theory fails qualitatively when N — K 
because the special points induce orbit closure in order 2^ steps. However, if local 
structures prevent the system from ever reaching a special point, then the mechanism 
for closing the orbits in order 2^ steps does not operate and it is plausible that the 
typical orbit lengths will be much longer than 2^. We test this idea by calculating 
orbit lengths in systems with one spin assigned the function —1 (so it follows the se- 
quence (+1,-1-1,— 1,-1, +1,-1-1, •••)) and a second individual spin assigned the func- 
tion +1 and the initial condition (+1, —1) (so its evolution is (+1, —1, +1, —1, ■■■))• 
Such (+, +, — , — ), (+, — , +, — ) systems can never reach a special point. As figure 8 
demonstrates, in these systems orbit lengths grow with N much faster than 2^; the 
numerical results are consistent with 2^^ scaling. 

We beheve that the orbit lengths in systems with N » 2^^ and randomly chosen 
functions cannot grow faster than 2'^^^^^^^^'^\ where e{K) is of order 2~^ . This is 
because of order N2~'^ spins have input functions 1 or —1 and hence cycle with 
periods 1, 2, and 4. More generally, we expect interesting crossover phenomena to 
occur when K is both large and of the order of log2{log2N). Even for K — ?>, simu- 
lating systems with large enough N to permit numerical exploration of these effects 
is computationally prohibitive. 
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Fig. 9. Survival functions of regular cycle lengths plotted against 1/2^ for various K values 
and for N = 10. The solid lines are the power law fits for the survival functions. By 
inspection one can see for K G [1.4, 1.7] the survival functions follow power law fairly well 
for three decades. 

4.-2 Average cycle length versus N for different values of K 

When K is large, the random hopping approximation works well and the cycle length 
distribution scales as uo = 2^ . Figure 5 demonstrates that the distribution of cycle 
lengths does not change dramatically while K decreases until K is quite small. There- 
fore, it is reasonable to expect the average cycle length to increase exponentially with 
N when K is large. On the other hand, when K — the average cycle length is 
bounded above by 4. For K = 1, our simulations and analytic arguments indicate 
that it scales as logA^; the results will be presented in the companion publication. 
For values of K in the range [1.4, 1.7], the survival functions S{1) decay roughly as 
a power law over three decades in cycle length as shown in figure 9. This result 
is consistent with the presence of a phase transition. The companion paper on the 
behavior of the Hamming distance presents more compelling evidence for a phase 
transition at X ~ 1.6. 
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5 Summciry 



This paper addresses the dynamics of a Boolean network model of N elements with K 
inputs with time-reversible dynamics. We present the general setup of the model and 
introduce the concept of special points and the distinction between special cycles and 
regular cycles. The relation between special points and the properties of the cycles 
is demonstrated. We show that the numbers of special points and special cycles for 
each realization are proportional to a; = 2^, where N is the number of variables in 
the system. 

We determine the probability distribution of cycle lengths as well as the survival 
functions. In limiting case K — 0, the cycle length is bounded above by 4 and the 
probability that a cycle length is 4 approaches 1 as N increases. For K — N, within a 
random hopping approximation we calculate the survival functions for regular cycles 
and for special cycles, which agree with simulational data extremely well. 

Finally, we present the simulational results for survival functions for intermediate 
K values. A population of anomalously long regular cycles scaling as 2^^ is found 
for small K values and explained based on the notions of special points and local 
structures. The correlation between typical cycle length and the K values of the 
networks is studied, and we find that the typical cycle length increases logarithmically 
with N when K < 1.4, exponentially when K > 1.7, and following a power law when 
K falls in between; these results are compatible with the presence of a phase transition 
for K somewhere in the range of [1.4, 1.7]. 
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Appendix A: Some statistical properties of sandwich points 

The text discusses twin and sandwich symmetry points which induce closure of orbits 
in the reversible Kauffman model. Each type arises when the substate is such 
that each of the functions of these inputs takes on a target value. For a twin point, 
cl^i = (Tf-, so the target function is different at every time step. For a sandwich point, 
for all t the target function is = 1 for every i. 

To calculate orbit lengths, we need to compute the probability that a symmetry point 
of cither type occurs at each time t. For a given realization of couplings, the number 
of substates for which the functions take on a particular value can vary. Because the 
target value for the twin point is different for different substates, whereas the target 
values for the sandwich points are the same at all times, the statistics of the two 
types of symmetry points are different. 

The process we consider is one in which couplings and an initial condition are chosen, 
and then the system evolves in time. We assume that this time evolution yields a 
random sampling of all possible spin configurations or substates. At each time t we 
examine whether a symmetry point of either type has been reached. Let m be the 
number of substates for which = 1 for all j (the criterion for a sandwich point), 
and rrit be the number of spin configurations for which each function takes on the 
target value for a twin point at time t. At a time t, the probability of being at a 
twin point is rrit/uj, whereas the probability of being at a sandwich point is m/uj. 
On average it takes many trials before a special point is reached; the probability of 
having observed a sandwich point after T trials ~ S^i^t/^ ~ (m)T/uj = T/uj, 
where () is the average over realizations, whereas the probability of having observed 
a sandwich symmetry point is mT/uj. 

We wish to calculate the probability that a randomly chosen realization has a given 
value of m. Now each output takes on a given value with probability 1/2, so on average 
the probability that N outputs all have given target values is (1/2)^, implying that 
the realization average (m) = 2"^. However, if one of the functions happens to be 
= —1, then clearly there are no sandwich points. We wish to find Pxim), the 
fraction of all possible reahzations of the couphngs for a given K that yield each 
value of m. 

First we consider K = N. This case is particularly relevant because when K is large, 
essentially all the orbits close because of the symmetry points. Here, the functions can 
be viewed as mapping a given input substate into a randomly chosen output substate. 
Since there are cu — 2^^ possible output substates, of which one is the target, each 
input configuration has a probability 1/uj of having its output be the target. Pk=n{0), 
the probability that no input configurations has as its output the target configuration, 
is (1 — 2^^)^ , which, as N ^ oo, approaches 1/e. The probability that exactly one 
of the 2^ different inputs yields the target output is (2^)(2-^)(l - 1/e. 
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Similarly. PK=N{m), the probability that exactly m of the 2^ different inputs yields 
the target output is 



PK=N{rn) = 




^ (AT^oo, m<2^). (.1) 

mle 

Thus, Pk=n{'iti) is a Poisson distribution. Since Pk=n{iti) falls off quickly as m gets 
large, clearly it is consistent to assume that m <^2^. 

However, when K is finite and N is large enough, we expect the behavior of PK{m) to 
differ quahtatively from the K — N result. We expect that almost all configurations 
will have m — for any finite K as N ^ oo. We have argued before that if in a 
realization a spin is assigned a function that equals to constant -1, the realization has 
no sandwich point. Also, if two spins are assigned functions and such that 

for all inputs, there can be no sandwich point for the realization. There are many other 
possible mechanisms that lead to m = 0. Clearly, the probabihty that a reahzation 
has at least one spin function of -1 bounds below the probability that it has no 
sandwich point. Among 2*^^ ^ possible functions that can be assigned to one spin, one 
is —1 for all inputs. Assuming that all functions are equally likely to be picked and 
that the function choices for different spins are independent, the probability that no 
spin is assigned the constant function -1 is 

(1 - «^ exp(-^) (Ar»2(^")). 

Thus, Pft:(m = 0) is bounded by the probability that the realization has at least one 
function that is —1; or Pft'(m = 0) > (1 — exp(— ^^^)). This result imphes that 
whenever K is finite, in a large enough system sandwich points cause orbit closure 
only in a vanishingly small fraction of realizations. However, when K is not small, 
realizations with sandwich points are rare only when N is enormous (when AT > 22''). 



Appendix B: Relation between special points and cycles 

Here we prove the results used in section 2 that 1) each special cycle contains exactly 
two special points, 2) that cycles with two special points of the same kind have even 
cycle lengths, and 3) cycles with different kinds of special points have odd lengths. 

To prove that each special cycle contains two and only two special points, we first 
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consider a cycle of even length 2n, 



^05 Si, . . . , T,2n-1- 

Suppose there is a twin point in the cycle. By relabeling the cycle, we can get 

by definition. Now = since the cycle is time reversible, so that Eq = S2„_i 

when we take t = n — 1. Thus we find another twin point in this even-length cycle. 
If there is a sandwich point at n, then 

and S„_i-t = S„,+i+t, thus 

Here we find another sandwich point at 0. Similarly, when the cycle is odd, we will 
find a twin point in the presence of a sandwich point, and vice versa. 

By now we have proven that if there is a special point in the cycle, then there has 
to be another. The statement that the cycle length being odd or even depends on 
whether the special points are of the same kind, is also clear from the above argument. 

To finish the proof, we need to demonstrate that no orbit can contain more than 
two special points. Assume that an orbit of length L with more than two special 
points exists. Choose the origin of time so that = (one does this by placing 
a sandwich substate at t = or placing twin substates at t = 1 and t = L), and let 
P be the smallest value for which Sp_„ = Sp+„ for all rzQ by assumption, P < L/2. 
Then we must have simultaneously ^l-j = and Sp_„ = Sp+„. Letting j = P — n 
yields Si_p+„ = Sp.^ = Sp+„. Letting q = P + n, we obtain SL_2p+g = Sq. Thus 
the orbit period is L — 2P, which is strictly less than L, so we have a contradiction. 

We can also show that these two special points must be different from one another. 
Suppose not. Consider an orbit of length L, and choose the origin of time so that 
Sq and Sp are the same special point, with, by assumption, P < L — 1. Applying 
the map yields Sj = 5(p+j) for any j, so the orbit repeats after P steps. But this 
contradicts the assumption that L is strictly greater than P. 



If the special point at P is a twin point, then 2P is odd, otherwise it is even. 
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